# Load libraries, install if needed
pkgs <- c("tidyverse", "readxl", "tidylog", "RCurl", "devtools", "patchwork",
"viridis", "RColorBrewer", "here", "ggnewscale", "ggh4x")
if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){
install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
}
invisible(lapply(pkgs, library, character.only = T))
# Source code for map plots
# You need: # devtools::install_github("seananderson/ggsidekick") # not on CRAN; library(ggsidekick)
devtools::source_url("https://raw.githubusercontent.com/maxlindmark/pred-prey-overlap/main/R/functions/map-plot.R")
options(ggplot2.continuous.colour = "viridis")
pal_temp <- brewer.pal(n = 10, name = "Paired")[c(2, 6)]
pal_oxy <- brewer.pal(n = 10, name = "Paired")[c(8, 4)]
# Set path
home <- here::here()Plot conditional effects and grid predictions from 02
Intro
Plot conditional effects from predictions in script “02-fit-sdms-conditional.qmd”
Load packages & source functions
Read predictions
preds <- read_csv(paste0(home, "/output/data_conditional_02_sdms.csv")) |>
filter(oxy > 3 & oxy < 9) |>
filter(temp > 2 & temp < 14) |>
separate(group, c("species", "life_stage"), sep = "_", remove = FALSE) |>
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage)) |>
group_by(group) |>
mutate(est_sc = exp(est)/max(exp(est))) # for plotting all together in heatmapRows: 15000 Columns: 15
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): group
dbl (14): temp, oxy, quarter, temp_sc, temp_sq, oxy_sc, year, year_f, quarte...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 5,700 rows (38%), 9,300 rows remaining
filter: removed 2,418 rows (26%), 6,882 rows remaining
mutate: changed 6,882 values (100%) of 'species' (0 new NA)
changed 6,882 values (100%) of 'life_stage' (0 new NA)
group_by: one grouping variable (group)
mutate (grouped): new variable 'est_sc' (double) with 3,621 unique values and 0% NAgrid_preds <- read_csv(paste0(home, "/output/data_pred_grids_02_sdms.csv")) |>
separate(group, c("species", "life_stage"), sep = "_", remove = FALSE) |>
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage))Rows: 5660568 Columns: 32
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (4): substrate, month_year, ices_rect, group
dbl (28): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_d...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
mutate: changed 5,660,568 values (100%) of 'species' (0 new NA)
changed 5,660,568 values (100%) of 'life_stage' (0 new NA)Calculate weighted quantiles
wq <- grid_preds |>
group_by(group, quarter) |>
summarise("Median_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.5),
"10_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.1),
"25_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.25),
"75_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.75),
"90_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.9),
"Env_oxy" = median(oxy),
"Median_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.5),
"10_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.1),
"25_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.25),
"75_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.75),
"90_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.9),
"Env_temp" = median(temp)) |>
ungroup() |>
separate(group, c("species", "life_stage"), sep = "_", remove = FALSE) |>
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage)) group_by: 2 grouping variables (group, quarter)summarise: now 12 rows and 14 columns, one group variable remaining (group)ungroup: no grouping variablesmutate: changed 12 values (100%) of 'species' (0 new NA) changed 12 values (100%) of 'life_stage' (0 new NA)wq_annual <- grid_preds |>
group_by(group, year, quarter) |>
# Could use reframe here...
summarise("Median_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.5),
"1st_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.25),
"3rd_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.75),
"Env_oxy" = median(oxy),
"Median_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.5),
"1st_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.25),
"3rd_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.75),
"Env_temp" = median(temp)) |>
ungroup() |>
separate(group, c("species", "life_stage"), sep = "_", remove = FALSE) |>
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage)) group_by: 3 grouping variables (group, year, quarter)summarise: now 348 rows and 11 columns, 2 group variables remaining (group, year)ungroup: no grouping variablesmutate: changed 348 values (100%) of 'species' (0 new NA) changed 348 values (100%) of 'life_stage' (0 new NA)Plot conditional effects
# Oxygen
preds_oxy <- preds |>
filter(temp %in% sort(unique(preds$temp))[0.5*length(unique(preds$temp))])filter (grouped): removed 6,696 rows (97%), 186 rows remainingwi_range_oxy <- preds_oxy |>
left_join(wq |>
filter(quarter == 4) |>
dplyr::select(group, `10_oxy`, `25_oxy`, `75_oxy`, `90_oxy`), by = "group")filter: removed 6 rows (50%), 6 rows remainingleft_join: added 4 columns (10_oxy, 25_oxy, 75_oxy, 90_oxy) > rows only in x 0 > rows only in y ( 0) > matched rows 186 > ===== > rows total 186# TODO: make predictions for the exact value??
wi_median_oxy <- preds_oxy |>
left_join(wq |>
filter(quarter == 4) |>
dplyr::select(group, Median_oxy), by = "group") |>
filter(oxy > 0.985*Median_oxy & oxy < 1.015*Median_oxy) |>
as.data.frame()filter: removed 6 rows (50%), 6 rows remainingleft_join: added one column (Median_oxy) > rows only in x 0 > rows only in y ( 0) > matched rows 186 > ===== > rows total 186filter (grouped): removed 180 rows (97%), 6 rows remainingp1 <- ggplot() +
ggh4x::facet_grid2(life_stage ~ species, scales = "free", independent = "y") +
# NOTE: 80% Confidence interval
# 25th and 75th percentile
geom_ribbon(data = wi_range_oxy |> filter(oxy > `25_oxy` & oxy < `75_oxy`), fill = "black",
aes(oxy, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)),
color = NA, alpha = 0.2) +
geom_line(data = wi_range_oxy |> filter(oxy > `25_oxy` & oxy < `75_oxy`), aes(oxy, exp(est)),
color = "black", alpha = 0.5, linewidth = 1.5) +
# 10th and 90th percentile
geom_ribbon(data = wi_range_oxy |> filter(oxy > `10_oxy` & oxy < `90_oxy`), fill = "black",
aes(oxy, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)),
color = NA, alpha = 0.2) +
geom_line(data = wi_range_oxy |> filter(oxy > `10_oxy` & oxy < `90_oxy`), aes(oxy, exp(est)),
color = "black", alpha = 0.5, linewidth = 1) +
# NOTE: 85% Confidence interval
geom_ribbon(data = preds_oxy, fill = "black",
aes(oxy, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)), color = NA, alpha = 0.2) +
geom_line(data = preds_oxy, aes(oxy, exp(est)), alpha = 0.5, linewidth = 0.5, color = "black") +
# Median!
geom_point(data = wi_median_oxy, aes(Median_oxy, exp(est)),
color = "black", alpha = 0.5, size = 2.5) +
labs(x = "Oxygen (ml/L)", y = "Biomass density (kg/km<sup>2</sup>)") +
theme(legend.position = "bottom",
legend.spacing.y = unit(0.1, 'cm'),
legend.key = element_rect(fill = "grey95", inherit.blank = TRUE),
axis.title.y = ggtext::element_markdown(),
plot.margin = unit(c(0,0,0,0), "cm")) +
NULLfilter (grouped): removed 141 rows (76%), 45 rows remainingfilter (grouped): removed 141 rows (76%), 45 rows remainingfilter (grouped): removed 101 rows (54%), 85 rows remaining
filter (grouped): removed 101 rows (54%), 85 rows remaining# Temperature
preds_temp <- preds |>
filter(oxy %in% sort(unique(preds$oxy))[0.5*length(unique(preds$oxy))])filter (grouped): removed 6,660 rows (97%), 222 rows remainingwi_range_temp <- preds_temp |>
left_join(wq |>
filter(quarter == 4) |>
dplyr::select(group, `10_temp`, `25_temp`, `75_temp`, `90_temp`), by = "group")filter: removed 6 rows (50%), 6 rows remainingleft_join: added 4 columns (10_temp, 25_temp, 75_temp, 90_temp) > rows only in x 0 > rows only in y ( 0) > matched rows 222 > ===== > rows total 222# TODO: make predictions for the exact value??
wi_median_temp <- preds_temp |>
left_join(wq |>
filter(quarter == 4) |>
dplyr::select(group, Median_temp), by = "group") |>
filter(temp > 0.985*Median_temp & temp < 1.015*Median_temp) |>
as.data.frame()filter: removed 6 rows (50%), 6 rows remainingleft_join: added one column (Median_temp) > rows only in x 0 > rows only in y ( 0) > matched rows 222 > ===== > rows total 222filter (grouped): removed 216 rows (97%), 6 rows remainingwi_median_temp temp oxy quarter temp_sc temp_sq oxy_sc year year_f quarter_f
1 7.698927 5.7852 1 0.6633954 0.4400935 -0.2013043 1999 1999 1
2 8.024739 5.7852 1 0.7803021 0.6088714 -0.2013043 1999 1999 1
3 10.305422 5.7852 1 1.5986492 2.5556793 -0.2013043 1999 1999 1
4 9.653798 5.7852 1 1.3648358 1.8627766 -0.2013043 1999 1999 1
5 8.350551 5.7852 1 0.8972089 0.8049837 -0.2013043 1999 1999 1
6 8.350551 5.7852 1 0.8972089 0.8049837 -0.2013043 1999 1999 1
sal_sc depth_sc depth_sq group species life_stage est
1 0 0 0 cod_adult Cod Adult 3.4416694
2 0 0 0 cod_juvenile Cod Juvenile 1.4838354
3 0 0 0 plaice_adult Plaice Adult -1.4838181
4 0 0 0 plaice_juvenile Plaice Juvenile -6.4750500
5 0 0 0 flounder_adult Flounder Adult 4.8314757
6 0 0 0 flounder_juvenile Flounder Juvenile -0.5487353
est_se est_sc Median_temp
1 0.6249370 0.9116380 7.602427
2 0.8294039 0.6392413 8.129976
3 0.8646872 0.7731394 10.387263
4 1.2343494 0.5440757 9.573396
5 0.3224810 0.4976172 8.387172
6 0.5137997 0.4212632 8.333579p2 <- ggplot() +
ggh4x::facet_grid2(life_stage ~ species, scales = "free", independent = "y") +
# NOTE: 80% Confidence interval
# 25th and 75th percentile
geom_ribbon(data = wi_range_temp |> filter(temp > `25_temp` & temp < `75_temp`), fill = "black",
aes(temp, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)),
color = NA, alpha = 0.2) +
geom_line(data = wi_range_temp |> filter(temp > `25_temp` & temp < `75_temp`), aes(temp, exp(est)),
color = "black", alpha = 0.5, linewidth = 1.5) +
# 10th and 90th percentile
geom_ribbon(data = wi_range_temp |> filter(temp > `10_temp` & temp < `90_temp`), fill = "black",
aes(temp, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)),
color = NA, alpha = 0.2) +
geom_line(data = wi_range_temp |> filter(temp > `10_temp` & temp < `90_temp`), aes(temp, exp(est)),
color = "black", alpha = 0.5, linewidth = 1) +
# NOTE: 85% Confidence interval
geom_ribbon(data = preds_temp, fill = "black",
aes(temp, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)),
color = NA, alpha = 0.2) +
geom_line(data = preds_temp, aes(temp, exp(est)), alpha = 0.5, linewidth = 0.5, color = "black") +
# Median!
geom_point(data = wi_median_temp, aes(Median_temp, exp(est)),
color = "black", alpha = 0.5, size = 2.5) +
labs(x = "Temperature (°C)", y = "Biomass density (kg/km<sup>2</sup>)") +
theme(legend.position = "bottom",
legend.spacing.y = unit(0.1, 'cm'),
legend.key = element_rect(fill = "grey95", inherit.blank = TRUE),
axis.title.y = ggtext::element_markdown(),
plot.margin = unit(c(0,0,0,0), "cm")) +
NULLfilter (grouped): removed 178 rows (80%), 44 rows remainingfilter (grouped): removed 178 rows (80%), 44 rows remainingfilter (grouped): removed 138 rows (62%), 84 rows remaining
filter (grouped): removed 138 rows (62%), 84 rows remaining(p1 / p2) + plot_annotation(tag_levels = "A")
ggsave(paste0(home, "/figures/oxy_temp_conditional.pdf"), width = 17, height = 17, units = "cm")Plot weighted quantiles over time
t <- wq_annual |>
pivot_longer(c("Median_oxy", `1st_oxy`, `3rd_oxy`, "Env_oxy", "Median_temp", `1st_temp`, `3rd_temp`, "Env_temp")) |>
mutate(var = ifelse(grepl("oxy", name), "Oxygen", "Temperature"),
type = ifelse(grepl("Env", name), "Environment", "Biomass-weighted")) |>
separate(group, c("species", "life_stage"), sep = "_", remove = FALSE) |>
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage),
group2 = paste(species, life_stage, sep = " ")) |>
filter(name %in% c("Median_oxy", "Median_temp", "Env_oxy", "Env_temp"))pivot_longer: reorganized (Median_oxy, 1st_oxy, 3rd_oxy, Env_oxy, Median_temp, …) into (name, value) [was 348x13, now 2784x7]mutate: new variable 'var' (character) with 2 unique values and 0% NA new variable 'type' (character) with 2 unique values and 0% NAmutate: changed 2,784 values (100%) of 'species' (0 new NA) changed 2,784 values (100%) of 'life_stage' (0 new NA) new variable 'group2' (character) with 6 unique values and 0% NAfilter: removed 1,392 rows (50%), 1,392 rows remainingt_env <- t |>
filter(type == "Environment" & group == "cod_adult") |>
mutate(group2 = "Environment")filter: removed 1,276 rows (92%), 116 rows remainingmutate: changed 116 values (100%) of 'group2' (0 new NA)t2 <- bind_rows(t_env, t |> filter(!type == "Environment"))filter: removed 696 rows (50%), 696 rows remaining# Reorder to have Environment at the end
t2$group2 <- factor(t2$group2, levels = rev(unique(t2$group2)))
pal <- brewer.pal(name = "Dark2", n = 8)[3:8]
t2 |>
mutate(quarter2 = ifelse(quarter == 1, "Quarter 1", "Quarter4")) |>
ggplot(aes(year, value, color = group2, linetype = group2)) +
geom_line(alpha = 0.8) +
#facet_grid(quarter ~ var) +
ggh4x::facet_grid2(var ~ quarter2, scales = "free") +
labs(y = "Environmental variable", x = "Year", linetype = "", color = "") +
guides(color = guide_legend(nrow = 3,override.aes = list(linetype = c(rep(1, 6), 2),
size = 0.5)),
linetype = "none") +
scale_x_continuous(breaks = c(seq(1993, 2021, by = 5))) +
scale_linetype_manual(values = c(rep(1, 6), 2)) +
scale_color_manual(values = c(rev(pal), "grey40")) +
theme(legend.position = c(0.24, 0.36),
legend.spacing.y = unit(-2, 'cm'),
legend.spacing.x = unit(0, 'cm'),
legend.background = element_rect(fill = NA))mutate: new variable 'quarter2' (character) with 2 unique values and 0% NA
ggsave(paste0(home, "/figures/wm_temp_oxy.pdf"), width = 17, height = 14, units = "cm")Extra plot for Gotje presentation
Correlation and slope plot
# Correlation plot
dcor <- t |> filter(!type == "Environment") |>
mutate(id = paste(name, year, quarter, sep = "_")) |>
dplyr::select(group2, year, var, value, quarter, id)filter: removed 696 rows (50%), 696 rows remainingmutate: new variable 'id' (character) with 116 unique values and 0% NAt_env2 <- t_env |>
mutate(name = ifelse(name == "Env_oxy", "Median_oxy", "Median_temp")) |>
mutate(id = paste(name, year, quarter, sep = "_")) |>
rename(env_value = value) |>
dplyr::select(id, env_value)mutate: changed 116 values (100%) of 'name' (0 new NA)
mutate: new variable 'id' (character) with 116 unique values and 0% NArename: renamed one variable (env_value)dcor2 <- dcor |> left_join(t_env2, by = "id")left_join: added one column (env_value) > rows only in x 0 > rows only in y ( 0) > matched rows 696 > ===== > rows total 696cors <- plyr::ddply(dcor2, c("group2", "var", "quarter"),
summarise, cor = round(cor(env_value, value), 2)) |>
arrange(var)summarise: now one row and one column, ungroupedsummarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungroupedhist(cors$cor)
# Plot correlations between environment and weighted
cor_pal <- brewer.pal(n = 8, name = "Dark2")[6:7]
dcor2 |>
mutate(quarter2 = ifelse(quarter == 1, "Quarter 1", "Quarter 4")) |>
ggplot(aes(env_value, value, color = quarter2)) +
ggh4x::facet_grid2(var ~ group2, scales = "free") +
geom_abline(color = "grey30", linetype = 2, linewidth = 0.5) +
geom_point(alpha = 0.7, size = 0.5) +
geom_text(data = cors |> filter(quarter == 1),
aes(label = paste("r=", cor, sep = "")), x = -Inf, y = Inf, hjust = -0.1, vjust = 1.5,
inherit.aes = FALSE, fontface = "italic", color = cor_pal[2], size = 2.5) +
geom_text(data = cors |> filter(quarter == 4),
aes(label = paste("r=", cor, sep = "")), x = -Inf, y = Inf, hjust = -0.1, vjust = 2.7,
inherit.aes = FALSE, fontface = "italic", color = cor_pal[1], size = 2.5) +
labs(y = "Biomass-weighted", x = "Environment", color = "") +
stat_smooth(method = "lm", se = FALSE, size = 0.5) +
scale_color_manual(values = cor_pal) +
scale_x_continuous(breaks = c(4:8)) +
theme_sleek(base_size = 10) +
theme(aspect.ratio = 1,
legend.key.size = unit(0.2, 'cm'),
legend.text = element_text(size = 6),
legend.position = c(0.11, 0),
legend.spacing.y = unit(-1, 'cm'))mutate: new variable 'quarter2' (character) with 2 unique values and 0% NAfilter: removed 12 rows (50%), 12 rows remaining
filter: removed 12 rows (50%), 12 rows remainingWarning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.`geom_smooth()` using formula = 'y ~ x'
ggsave(paste0(home, "/figures/supp/wm_cors.pdf"), width = 17, height = 7, units = "cm")`geom_smooth()` using formula = 'y ~ x'# Plot the slopes over time for group and quarter, compare to the environment slope. Boxplots and lines??
s_slopes <- t |> filter(!type == "Environment") |>
mutate(id = paste(group, var, quarter, sep = "_")) |>
ungroup() |>
dplyr::select(year, id, value)filter: removed 696 rows (50%), 696 rows remainingmutate: new variable 'id' (character) with 24 unique values and 0% NAungroup: no grouping variablest_slopes <- t_env |>
mutate(id = paste("Environment", group2, var, quarter, sep = "_")) |>
dplyr::select(year, id, value)mutate: new variable 'id' (character) with 4 unique values and 0% NAslopes <- bind_rows(s_slopes, t_slopes)
# Calculate slopes over time
slopes2 <- slopes %>%
split(.$id) %>%
purrr::map(~lm(value ~ year, data = .x)) |>
purrr::map_df(broom::tidy, .id = 'Year') |>
filter(term == "year") |>
mutate(upr = estimate + 1.96*std.error,
lwr = estimate - 1.96*std.error,
sig = ifelse(estimate > lwr & estimate < upr, "sig.", "not sig.")) |>
rename(id = Year,
year_slope = estimate) |>
dplyr::select(id, year_slope, sig, lwr, upr) |>
separate(id, into = c("species", "life_stage", "variable", "quarter"), remove = FALSE) |>
#mutate(group3 = paste(str_to_title(species), str_to_title(life_stage), paste0("Q", quarter)))
mutate(group3 = paste(str_to_title(species), str_to_title(life_stage)),
x = "x") |>
mutate(quarter2 = ifelse(quarter == 1, "Quarter 1", "Quarter 4"))filter: removed 28 rows (50%), 28 rows remainingmutate: new variable 'upr' (double) with 28 unique values and 0% NA new variable 'lwr' (double) with 28 unique values and 0% NA new variable 'sig' (character) with one unique value and 0% NArename: renamed 2 variables (id, year_slope)mutate: new variable 'group3' (character) with 7 unique values and 0% NA new variable 'x' (character) with one unique value and 0% NAmutate: new variable 'quarter2' (character) with 2 unique values and 0% NAhlines <- slopes2 |>
filter(species == "Environment")filter: removed 24 rows (86%), 4 rows remainingpal <- brewer.pal(name = "Dark2", n = 8)[3:8]
slopes2 |>
filter(!species == "Environment") |>
ggplot(aes(x, year_slope, ymin = lwr, ymax = upr, color = group3)) +
geom_rect(data = hlines, aes(xmin = -Inf, xmax = Inf, ymin = lwr, ymax = upr, fill = "Env. slope 95% CI"),
inherit.aes = FALSE, alpha = 0.1) +
geom_hline(yintercept = 0, linetype = 3, color = "tomato2", alpha = 1) +
geom_hline(data = hlines, aes(yintercept = year_slope, linetype = "Env. slope"), alpha = 0.7) +
geom_point(size = 2.5, position = position_dodge(width = 0.4)) +
geom_errorbar(position = position_dodge(width = 0.4), width = 0) +
ggh4x::facet_grid2(variable ~ quarter2, scales = "free") +
scale_fill_manual(values = "grey10") +
scale_linetype_manual(values = 2) +
scale_color_manual(values = pal) +
labs(color = "", y = "Year-slope", linetype = "", fill = "") +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
legend.position = "bottom") +
NULLfilter: removed 4 rows (14%), 24 rows remaining
ggsave(paste0(home, "/figures/year_slopes.pdf"), width = 17, height = 17, units = "cm")Plot spatially-varying quarter effects, spatial predictions, and ST random fields
grid_preds <- grid_preds |>
mutate(group2 = paste(species, life_stage, sep = " ")) mutate: new variable 'group2' (character) with 6 unique values and 0% NAnames(grid_preds) [1] "X" "Y" "year"
[4] "lon" "lat" "depth"
[7] "substrate" "quarter" "month"
[10] "month_year" "oxy" "temp"
[13] "sal" "sub_div" "ices_rect"
[16] "temp_sc" "temp_sq" "oxy_sc"
[19] "oxy_sq" "sal_sc" "depth_sc"
[22] "depth_sq" "quarter_f" "year_f"
[25] "est" "est_non_rf" "est_rf"
[28] "omega_s" "zeta_s_quarter_f1" "zeta_s_quarter_f4"
[31] "epsilon_st" "group" "species"
[34] "life_stage" "group2" # Plot spatially-varying effect
sv <- grid_preds |>
#filter(group == pull(grid_preds, group)[1]) |>
filter(year == 1999) |>
pivot_longer(c("zeta_s_quarter_f1", "zeta_s_quarter_f4"),
names_to = "field", values_to = "zeta") |>
mutate(Quarter = ifelse(field == "zeta_s_quarter_f1", "Quarter 1", "Quarter 4"))filter: removed 5,465,376 rows (97%), 195,192 rows remainingpivot_longer: reorganized (zeta_s_quarter_f1, zeta_s_quarter_f4) into (field, zeta) [was 195192x35, now 390384x35]mutate: new variable 'Quarter' (character) with 2 unique values and 0% NAplot_map_fc +
geom_raster(data = sv, aes(X*1000, Y*1000, fill = zeta)) +
facet_grid(Quarter ~ group2) +
scale_fill_gradient2(name = "Zeta")Warning: Removed 18264 rows containing missing values (`geom_raster()`).
ggsave(paste0(home, "/figures/supp/sv_effects.pdf"), width = 19, height = 9, units = "cm")Warning: Removed 18264 rows containing missing values (`geom_raster()`).for(i in unique(grid_preds$group2)){
dd <- grid_preds |> filter(group2 == i)
j <- pull(dd, group)[1]
# Quarter 1
print(
plot_map_fc +
geom_raster(data = filter(dd, quarter == 1), aes(X*1000, Y*1000, fill = exp(est))) +
facet_wrap(~year) +
scale_fill_viridis(trans = "sqrt", name = "Biomass density (kg/km)",
# trim extreme high values to make spatial variation more visible
na.value = "yellow", limits = c(0, exp(quantile(filter(dd, quarter == 4)$est, 0.999)))) +
labs(caption = paste("maximum estimated biomass density =", round(max(exp(filter(dd, quarter == 4)$est)))),
title = i, subtitle = "Quarter 1")
)
ggsave(paste0(home, paste0("/figures/supp/spatial_biomass", "_q1_", j, ".pdf")), width = 17, height = 17, units = "cm")
# Quarter 4
print(
plot_map_fc +
geom_raster(data = filter(dd, quarter == 4), aes(X*1000, Y*1000, fill = exp(est))) +
facet_wrap(~year) +
scale_fill_viridis(trans = "sqrt", name = "Biomass density (kg/km)",
# trim extreme high values to make spatial variation more visible
na.value = "yellow", limits = c(0, exp(quantile(filter(dd, quarter == 4)$est, 0.999)))) +
labs(caption = paste("maximum estimated biomass density =", round(max(exp(filter(dd, quarter == 4)$est)))),
title = i, subtitle = "Quarter 4")
)
ggsave(paste0(home, paste0("/figures/supp/spatial_biomass", "_q4_", j, ".pdf")), width = 17, height = 17, units = "cm")
}filter: removed 4,717,140 rows (83%), 943,428 rows remainingfilter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remainingWarning: Removed 22069 rows containing missing values (`geom_raster()`).Warning: Removed 22069 rows containing missing values (`geom_raster()`).filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).filter: removed 4,717,140 rows (83%), 943,428 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).filter: removed 4,717,140 rows (83%), 943,428 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).filter: removed 4,717,140 rows (83%), 943,428 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).filter: removed 4,717,140 rows (83%), 943,428 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).filter: removed 4,717,140 rows (83%), 943,428 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).